Lesson 6 - Scaling up your Analyses: Writing Functions in R

---

A quick intro to the intro to R Lesson Series

This ‘Intro to R Lesson Series’ is brought to you by the Centre for the Analysis of Genome Evolution & Function’s (CAGEF) bioinformatics training initiative. This course was developed based on feedback on the needs and interests of the Department of Cell & Systems Biology and the Department of Ecology and Evolutionary Biology.

This lesson is the last in a 6-part series. Hopefully you have learned a bit about how to import and manipulate data, make exploratory plots, perform some basic statistical tests, test a regression model, and make some even prettier plots and documents to share to results.

Today we are going to learn to write functions, which can save you time and help scale up your analyses. It will also help to minimize copy and pasting, redundancy, and feeling like a newb.

The structure of the class is a code-along style. It is hands on. The lecture AND code we are going through are available on GitHub for download at https://github.com/eacton/CAGEF, so you can spend the time coding and not taking notes. As we go along, there will be some challenge questions and multiple choice questions on Socrative. At the end of the class if you could please fill out a post-lesson survey (https://www.surveymonkey.com/r/X5C2JGR), it will help me further develop this course and would be greatly appreciated.

---

Packages Used in This Lesson

The following packages are used in this lesson:

tidyverse (ggplot2, tidyr, dplyr)
testthat assertthat

Please install and load these packages for the lesson. In this document I will load each package separately, but I will not be reminding you to install the package. Remember: these packages may be from CRAN OR Bioconductor.

---

Highlighting

grey background - a package, function, code or command
italics - an important term or concept
bold - heading or a term that is being defined
blue text - named or unnamed hyperlink

---

Objective: At the end of this session you will be able to write functions, making your coding more efficient, understandable, reproducible, and hopefully less frustrating.

Load the packages!

library(tidyverse)
library(testthat)
library(assertr)

Error in library(assertr) : there is no package called ‘assertr’

Dataset: Pyrosequencing of the V3-V5 hypervariable regions of the 16S rRNA gene

16S rRNA pyrosequencing of 30 latrines from Tanzania and Vietnam at different depths (multiples of 20cm). Microbial abundance is represented in Operational Taxonomic Units (OTUs). Operational Taxonomic Units (OTUs) are groups of organisms defined by a specified level of DNA sequence similarity at a marker gene (e.g. 97% similarity at the V4 hypervariable region of the 16S rRNA gene). Intrinsic environmental factors such as pH, temperature, organic matter composition were also recorded.

We have 2 csv files:

1. A metadata file (Naming conventions: [Country_LatrineNo_Depth]) with sample names and environmental variables.
   
2. OTU abundance table.

B Torondel, JHJ Ensink, O Gunvirusdu, UZ Ijaz, J Parkhill, F Abdelahi, V-A Nguyen, S Sudgen, W Gibson, AW Walker, and C Quince. Assessment of the influence of intrinsic environmental and geographical factors on the bacterial ecology of pit latrines Microbial Biotechnology, 9(2):209-223, 2016. DOI:10.1111/1751-7915.12334

---

Why Functions?

dat <- read.csv("data/long_SPE_pitlatrine.csv", header = T, stringsAsFactors = F)

I want to know the median and median absolute deviation of OTUs in all of our samples for Clostrida.

Great. What about the median and median absolute deviation for Bacilli? Hmmmm. I guess I’ll just copy and paste the above code because it will take 2 seconds.

Maybe I’ll check out Bacteroidia, too… by copy and pasting AGAIN.

Come to think of it, there are a lot of zeros in this data. Maybe I should be looking at mean and standard deviation instead of median and mad. I’ll just go back and change this in each of my 3 lines of code.

dat %>% filter(Taxa == "Clostridia") %>% summarize(mean = mean(OTUs), sd = sd(OTUs))
dat %>% filter(Taxa == "Bacilli") %>% summarize(mean = mean(OTUs), sd = sd(OTUs))
dat %>% filter(Taxa == "Bacteroidia") %>% summarize(mean = mean(OTUs), sd = sd(OTUs))

Does any of this feel familiar??

Functions allow us to:

• avoid repetition, make our code more compact
• avoid errors, copy/paste errors, typos
• versatile, updateable
• upscale, generalizable
• resuable, even for other projects!

The Structure of a Function

A function has 3 parts:

• formals are the input(s) of the function,
• the body is the braced expression,
• and the environment is the enclosure in which the function was defined.

We can look at these for a function we are familiar with, read.csv. Each argument is specified with $ and any default values are listed [1]. To use read.csv you must input a file - this argument has no default value. The rest of the arguments have a default value (such as header = TRUE), but can be modified by specifying the argument (by inputting header = FALSE). Most function defaults as ‘sensible’, but if you have not used the function before it is worth reading about the arguments in the help menu. ESPECIALLY do this with mathematical functions to see what assumptions you are making by using the default parameters. read.csv also has the magical ... argument. This argument allows you to pass arbitrary arguments to another function. We will see this in use shortly.

formals(read.csv)

$file


$header
[1] TRUE

$sep
[1] ","

$quote
[1] "\""

$dec
[1] "."

$fill
[1] TRUE

$comment.char
[1] ""

$...

The body of read.csv looks suspiciously like the function read.table because it is, indeed a ‘wrapper’ function which uses an existing function with different values for default arguments.

body(read.csv)

read.table(file = file, header = header, sep = sep, quote = quote, 
    dec = dec, fill = fill, comment.char = comment.char, ...)

The environment for read.csv is the utils package. When we make a function, it will be created in ‘.GlobalEnv’ the ‘global environment’. The global environment is our interactive workspace, viewed in the environment window in the top right, and can be thought of as a bag that holds all of our variables and functions. Note that we can use the dropdown menu to change the environment to package:utils, and we can find read.csv in that environment.

environment(read.csv)

<environment: namespace:utils>

Writing Functions

The basic form of a function takes in arguments (specifications by the user), has statements or expressions (what the function does), and returns an object. A function does not need to be named to be used, however, the greatest utility comes from reusing functions. A named function is called using parentheses: named_function(arguments).

Different from other programming languages, R is not picky about whitespace (the arrangement of tabs, spaces, indents). However formatting your function as below (statements indented below your arguments) will help to make it readable and familiar to other languages.

This may be a good time to point out that placing your cursor beside a closing bracket will highlight its corresponding opening bracket. This is really useful as your code gets more complex. Also, try deleting the closing bracket beside ‘object’. You will see that R is giving you a series of warnings (brackets highlighted in red, and ‘x’ beside 3 lines of code). Again, this is a helpful warning that saves you time from running an incomplete function.

named_function <- 
  function(arguments) {
      statments
      return(object) 
  }

We want to write a function that will calculate the mean and standard deviation for a given Taxa.

We can start by naming our function something vaguely informative. (Technically, you can name your functions and variables whatever you want. I could name them after Avengers characters, but people (or myself in 3 months) would have to work pretty hard to figure out what I was doing. Remember, with great power comes great responsibility).

Our function is going to take ‘Taxa’ as an argument. For some of you, when you type ‘{’ R will automatically give you your closing bracket ‘}’ so you won’t forget - ever trying to be helpful, thanks R. Next we can take our previous code, and save the output into an object called ‘summary_dat’. This is because I don’t want to overwrite my original data frame every time I run the function. We replace our specific Taxa, ‘Clostridia’, with our argument, ‘Taxa’. This is our expression. All we have to do now is return our data frame and make sure all of our brackets match.

Voila!

mean_func <- 
  function(Taxa) {
    summary_dat <- filter(dat, Taxa == Taxa) %>% 
    summarize(mean = mean(OTUs), sd = sd(OTUs))
    return(summary_dat)
}

When you run this function, it will appear, named and with its arguments, in the global environment. If you click on the white box to the right of the function, you can see the code for your function. You can now use your function.

out <- mean_func("Bacteriodia")

[1] "This is the mean and sd for Bacteriodia"

mean_func <- 
  function(Taxa) {
    filter(dat, Taxa == Taxa) %>% 
    summarize(mean = mean(OTUs), sd = sd(OTUs))
    print(paste('This is the mean for', Taxa))
}

Understanding Lexical Scoping

• why does dat get changed, but summary_dat does not end up being a global environmental variable when returned?

-may want to be able to take into different data frames…

mean_func <- function(dat, x) {
  dat <- filter(dat, Taxa == x) %>% 
    summarize(mean = mean(OTUs), sd = sd(OTUs))
}

mean_func(dat, "Clostridia")

Allows us to vary the Taxa and the functions.

mean_func("Clostridia", funs = list(mean = mean, median = median, sd = mean))

$mean
[1] 3423.407

$median
[1] 1432

$sd
[1] 3423.407

Defensive Programming

Testing Arguments and Validity of your Function

You might think that you are done with this function. You have the answers, right? But functions are not foolproof.

mean_func("Clostridia", funs = list(mean = quantile, median = median, sd = mean))

$mean
   0%   25%   50%   75%  100% 
    0   496  1432  4761 17572 

$median
[1] 1432

$sd
[1] 3423.407

Hmmm… mislabelled numbers could get a little confusing.

What if I wanted the summary statistics for 3 Taxa at once?

Uh oh. That’s bad. I didn’t want to merge the 3 Taxa… or did it take the last value?

Your next step is to try and intentionally break your function. Give it input that wasn’t meant for it and see what happens. Trouble-shooting this way is called ‘defensive programming’. What we are going to do is implement checks and balances such that our function can’t do anything weird.

Initial checks can be things like - what happens when I give an input that is numeric when I am expecting character data? What happens if there is an NA in the dataset? What happens if there is a typo? What happens if I am given an unexpected data structure? What happens if there is more or less data than I am expecting? What happens when I input logical data?

mean_func()

Error Handling

Unit Tests

-testthat assumption for function name = name?

Returns the last line… unless explicitly called with return

test_that("input is Taxa", {
  expect_success(expect_type(mean_func("Bacilli"), type = "character"))
  expect_failure(expect_type(mean_func(1), type = "character"))
})

[1] "This is the mean for Bacilli"
[1] "This is the mean for 1"

Error: Test failed: 'input is Taxa'
* Expectation did not fail

Assertions

In a perfect world, if something goes wrong while you are coding, you want an informative error message to reduce time trouble-shooting. We also want to do that for other people. Again, Hadley has come to the rescue with a way to test assertions and give meaningful error messages with the package assertthat.

Assertions are statements that are known to be TRUE, but we use whether these statements are TRUE or FALSE to our advantage in testing. For example, you can test whether the input to a function is a specific data type, length, or other custom argument by placing the assertion inside the function. You can also customize the output messages.

Let’s write some assertions into our function to test that:

-our input is a character vector -our input is in our list of Taxa (no typos!) -our input is only 1 value

We will then see how well this takes care of all of our inputs that ‘broke’ the function earlier. We will assess the error messages to see if they are clear.

function(Taxa) {
    assert_that(is.character(Taxa)) %>%
    validate_that(Taxa %in% dat$Taxa) %>%
    assert_that(length(Taxa) == 1) %>%
    summary <- filter(dat, Taxa == Taxa) %>% 
    summarize(mean = mean(OTUs), sd = sd(OTUs))
    return(summary)
    print(paste('This is the mean for', Taxa,'. ', summary$mean))
    
}

function(Taxa) {
    assert_that(is.character(Taxa)) %>%
    validate_that(Taxa %in% dat$Taxa) %>%
    assert_that(length(Taxa) == 1) %>%
    summary <- filter(dat, Taxa == Taxa) %>% 
    summarize(mean = mean(OTUs), sd = sd(OTUs))
    return(summary)
    print(paste('This is the mean for', Taxa,'. ', summary$mean))
    

}

I usually test to make sure it still works on something I expect it to when I make changes to my code.

mean_func(c("Clostridia", "Bacteriodia", "Bacilli"))

Error: assert_that: length of assertion is not 1

My function is still working on expected input.

For a vector of length 3, we receive the error that the ‘length of the assertion is not 1’. This could be more explicit for a user ie. “Please input only one Taxa”, but it is not bad. Note that using assert_that() we are not getting an output. The function stops when it gets a value that is FALSE and does not continue evaluation. This is much better than previously when we still got an output that was incorrect and did not get any warning. This also means that this passed the assertions for being a character vector and being in the Taxa in our dataset without error.

mean_func("Bacilllli")

Error: Bacilllli is not a Taxa in the dataset.

mean_func()

Error in eval(assertion, env) : 
  argument "Taxa" is missing, with no default

Most of these error messages are informative, but maybe not as explicit as they could be. Let’s customize our messages.

mean_func <- function(Taxa) {
    #assertions to test input
    assert_that(is.character(Taxa), msg = paste("Taxa is expected to be of type character. You have input type", typeof(Taxa), ".")) 
    assert_that(Taxa %in% dat$Taxa, msg = paste(Taxa, "is not a Taxa in the dataset."))
    assert_that(length(Taxa) == 1, msg = "Please input only one Taxa.")
    summary <- filter(dat, Taxa == Taxa) %>% 
    summarize(mean = mean(OTUs), sd = sd(OTUs))
 
    print(paste('The mean for', Taxa,'is', round(summary$mean,2), "."))
    return(summary)
    

}

We could be even more helpful and show the option of all Taxa to choose from, or convert everything to lowercase to reduce the chance of errors from having a mix of uppercase and lowercase. These are things you should think about, especially if you are passing on your code or making a package yourself. The easier people can use your functions, the more adoptable they will be.

There are a few functions included in assertthat for commonly used assertions. is.string() looks for a character vector of length 1 (combining 2 of our above assertions). Similarily is.count() looks for a single positive integer. noNA() tests the assertion that there are no NAs (different from is.na() because assertions are TRUE). Remember that the order of your assertions will affect the error message you receive. For example, an NA is not a string, but this fact may be masked if the check for a string precedes the check for an NA since the function will stop at the first FALSE evaluation.

not_empty() tests that all dimensions are greater than 0.

Let’s change a few assertions to test for an empty function call, NA as an input, and the use of is.string().

mean_func <- function(Taxa) {
    
    assert_that(not_empty(Taxa))
    assert_that(noNA(Taxa)) 
    assert_that(is.string(Taxa)) 
    assert_that(Taxa %in% dat$Taxa, msg = paste(Taxa, "is not a Taxa in the dataset.")) 
    
    summary <- filter(dat, Taxa == Taxa) %>% 
    summarize(mean = mean(OTUs), sd = sd(OTUs))
 
    print(paste('The mean for', Taxa,'is', round(summary$mean,2)))
    return(summary)
    

}

mean_func()

Error in dim(x) %||% length(x) : 
  argument "Taxa" is missing, with no default

We are told there is a missing value, but not where is is. (case for data frames?)

There are more of these helper functions to assert_that() which I encourage you to check out in the help menu.

Multiple Inputs

So far, we have seen functions that can be reused for different inputs. What if we want to run the function for many inputs? What if I want the function to be able to take 20 Taxa? We need to be able to iterate over multiple inputs, one at a time, in a controlled manner. Enter control statements. With control statements, an expression is set that if evaluated to be TRUE will continue proceed and if evaluated to be FALSE will not.

For loops will iterate through a specified number of items. Starting with a simple loop to print the numbers 1 through 10. This for statement will iterate through the vector of numbers 1:10. In this case i will be 1,2,3,4,5,6,7,8,9, then 10.

for (i in 1:10){
  print(i)
}

[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6
[1] 7
[1] 8
[1] 9
[1] 10

Often, we want to iterate over the contents of a vector, rather than a number. This is really the same thing.

for (i in x){
  print(i)
}

[1] "rabbit"
[1] "dog"
[1] "llama"

You can subset for the actual value.

for (i in 1:length(x)){
  print(x[i])
}

[1] "rabbit"
[1] "dog"
[1] "llama"

Or use numbers.

for (i in 1:length(x)){
  print(i)
}

[1] 1
[1] 2
[1] 3

Why use numbers?

for (i in 1:length(x)){
  if (i < 3){
  print(i)
  }
}

[1] 1
[1] 2

for (i in dat$Taxa){
  print(i)
  break
}

[1] "Acidobacteria_Gp1"

for (i in 1:nrow(dat)){
  print(dat[i,1])
  break
}

[1] "Acidobacteria_Gp1"

mean_func(c("Clostridia", "Bacteroidia", "Bacilli"))

[1] "The mean for Clostridia is 3423.41"
[1] "The mean for Bacteroidia is 684.57"
[1] "The mean for Bacilli is 330.85"

mean_func <- function(x, funs = list(mean = mean, sd = sd)) {
  dat <- filter(dat, Taxa == x) 
  lapply(funs, function(f) f(dat$OTUs))
  return(dat)
}

load_packages(c("foreach","MASS","doParallel"))

Loading required package: iterators
Loading required package: parallel

generalize to country, bacteria, #otus above/below mean? ranking position?

Basic Syntax of Functions in R

• the structure of a function (formals, body, parent environment)
• understanding scoping
• naming conventions
• generalizing a function
• setting default values
• adding ‘…’
• the output of your function
• sourcing functions
• NAs in functions
• return

Control Structures - if else - for - while - repeat - break - next

Handling Errors

• stopifnot(), warning(), suppressWarnings(), tryCatch()
• if stop

Testing Arguments and Validity of your Function

• informal testing
• testthat() for formal unit testing
• assertthat

Upscaling with functions

functions can be:

• assigned in variables
• stored in lists
• passed as arguments to other functions
• created inside other functions
• output by other functions

DRY principle - Don’t Repeat Yourself

---

Challenge

Write a function to check that the packages c(“dplyr”, “readxml”, “tidyr”) are installed and to load the packages. The function should install the packages if they are not already installled. Write a warning if the package is not installed, but is being installed. Will this function work for other packages?

---

---

Challenge

Write a function to read all .csv files in a directory into R and save each of them in a separate data frame.

read_files(path, FALSE)

[[1]]

[[2]]

[[3]]

NA

---

Resources:

http://stat545.com/block011_write-your-own-function-01.html
http://stat545.com/block011_write-your-own-function-02.html
http://stat545.com/block011_write-your-own-function-03.html
http://mazamascience.com/WorkingWithData/?p=912 http://adv-r.had.co.nz/Functions.html

Post-Lesson Assessment

---

Your feedback is essential to help the next cohort of trainees. Please take a minute to complete the following short survey: https://www.surveymonkey.com/r/X5C2JGR

---

Thanks for coming!!!